%% Exercise on Symmetrical Components

clearvars
close all;
clc;

% For the medium-voltage network shown in Figure 1, determine the voltages in the phasor domain at point 1 and point 2 using the theory of symmetrical components.
%
% Figure 1 - Single line diagram of the study network.
% Given data
% Nominal voltage of the source

Vn=15000;   %[V] RMS L-L

% Source short-circuit power

Scc=1.5E9;  %[VA]

% Relation between the real and imaginary parts of the short circuit impedance of the source \frac{R_{\textrm{cc}} }{X_{\textrm{cc}} }=0\ldotp 1.
% The power supply system is symmetrical three-phase with only the direct component (zero homopolar and inverse components). The three simple voltages can therefore be expressed in terms of the nominal compound voltage.

Ea_mod=Vn/sqrt(3)*sqrt(2);  %[V] amplitude L-G
Eb_mod=Vn/sqrt(3)*sqrt(2);  %[V]
Ec_mod=Vn/sqrt(3)*sqrt(2);  %[V]
Ea_ph=0;            %[rad]
Eb_ph=2/3*2*pi;       %[rad]
Ec_ph=1/3*2*pi;       %[rad]

Ea=Ea_mod*exp(1j*Ea_ph);
Eb=Eb_mod*exp(1j*Eb_ph);
Ec=Ec_mod*exp(1j*Ec_ph);

radtodeg=180/pi;    % conversion radians to degrees
degtorad=pi/180;    % conversion degrees to radians

spplot([Ea_mod Ea_ph*radtodeg],[Eb_mod Eb_ph*radtodeg],[Ec_mod Ec_ph*radtodeg])
legend('','','','','Ea','Eb','Ec');

% The load current system is unbalanced and asymmetrical:

Ia_mod=120;  %[A]
Ib_mod=90;   %[A]
Ic_mod=150;  %[A]

Ia_ph=27*degtorad;        %[rad]
Ib_ph=225*degtorad;       %[rad]
Ic_ph=125*degtorad;       %[rad]

Ia=Ia_mod*exp(1j*Ia_ph);
Ib=Ib_mod*exp(1j*Ib_ph);
Ic=Ic_mod*exp(1j*Ic_ph);

spplot([Ia_mod Ia_ph*radtodeg],[Ib_mod Ib_ph*radtodeg],[Ic_mod Ic_ph*radtodeg])
legend('','','','','Ia','Ib','Ic');

% Length of the line

L=5;     %[km]

% Cross-section of the line

xs_line=25;                %[mm^2]

% Resistance per unit length r_{l,\textrm{pul}}, reactance per unit length x_{l,\textrm{pul}} and impedance \bar{Z_l } =L\left(r_{l,\textrm{pul}} +j*x_{l,\textrm{pul}} \right) of the line

Rl_pu=0.929;             %[Ohm/km]
Xl_pu=0.15;              %[Ohm/km]
Zl=(Rl_pu+1i*Xl_pu)*L    %[Ohm]
Zl_mod=abs(Zl)           %[Ohm]
Zl_ph=angle(Zl)          %[rad]

%
% Question
% Determine the voltages in the phasor domain for each phase in buses #1 and #2.
%
% Solution
% Since the load is unbalanced and asymmetrical, in order to solve the three-phase circuit, the original circuit must be decomposed using symmetrical components.
%
% Figure 2. Circuits for each of the three sequences
% Short-circuit impedance of the source Z_{\textrm{cc}} =\frac{V_n^2 }{S_{\textrm{cc}} }

Zcc_mod=Vn^2/Scc   %[Ohm]

% Knowing the ratio between the real and imaginary parts of the source's short-circuit impedance, we can write the following relationships: \left(Z_{\textrm{cc}}^2 \right)=\left(R_{\textrm{cc}}^2 \right)+\left(X_{\textrm{cc}}^2 \right)=1\ldotp 01X_{\textrm{cc}}^2
% X_{\textrm{cc}} =\frac{Z_{\textrm{cc}} }{\sqrt{1\ldotp 01}}

Xcc=Zcc_mod/sqrt(1.01)      %[Ohm]
Rcc=Xcc*0.1                 %[Ohm]
Zcc=Rcc+1i*Xcc              %[Ohm]
Zcc_ph=phase(Zcc)           %[rad]

% The source and line impedance matrices are diagonal matrices equal to: {\bar{Z} }_{\textrm{cc}} I et {\bar{Z} }_l I
% where I is the identity matrix.

Zcc_mat=Zcc*eye(3)

Zl_mat=Zl*eye(3)

% Through the definition of the transformation matrix T_s =\left\lbrack \begin{array}{ccc}
% 1 & 1 & 1\\
% 1 & \alpha^2  & \alpha \\
% 1 & \alpha  & \alpha^2
% \end{array}\right\rbrack, it is easy to verify that the products {\left\lbrack T_s \right\rbrack }^{-1} \left\lbrack {\bar{Z} }_{\textrm{cc}} \right\rbrack \left\lbrack T_s \right\rbrack and {\left\lbrack T_s \right\rbrack }^{-1} \left\lbrack {\bar{Z} }_l \right\rbrack \left\lbrack T_s \right\rbrack are diagonal matrices whose diagonal components are the homopolar, direct and inverse sequence impedances.

alpha=exp(1j*2/3*pi)    % definition of alpha from the theory
Ts=[1 1 1; 1 alpha^2 alpha; 1 alpha alpha^2]
Zs_0di=Ts\Zcc_mat*Ts  % matrix multiplication
Zs_0di2 = Ts\Zcc_mat;
Zs_0=Zs_0di(1,1)
Zs_d=Zs_0di(2,2)
Zs_i=Zs_0di(3,3)

Zl_0di=Ts\Zl_mat*Ts
Zl_0=Zl_0di(1,1)
Zl_d=Zl_0di(2,2)
Zl_i=Zl_0di(3,3)

% Therefore {\bar{Z} }_0^S ={\bar{Z} }_i^S ={\bar{Z} }_d^S ={\bar{Z} }_{\textrm{cc}} and {\bar{Z} }_0^l ={\bar{Z} }_i^l ={\bar{Z} }_d^l ={\bar{Z} }_l (this result is expected because the two matrices \left\lbrack {\bar{Z} }_{\textrm{cc}} \right\rbrackand \left\lbrack {\bar{Z} }_l \right\rbrack are diagonal and with the same diagonal elements).
% Now we need to calculate the symmetrical components of the supply voltage and load currents. As long as the source is symmetrical and balanced, there is only the direct component in the triplet of supply voltages:
% {\bar{E} }_d^S =\frac{V_n }{\sqrt{3}}e^{\mathrm{j0}} ;{\bar{E} }_i^S ={\bar{E} }_0^S =0\ldotp

E_phases=[Ea; Eb; Ec]
E_seq=Ts\E_phases

Ed_mod=abs(E_seq(2))
Ed_ph=angle(E_seq(2))
Ed=E_seq(2)

% For the load currents we have:

I_phases=[Ia; Ib; Ic]
I_seq=Ts\I_phases

I0=I_seq(1)
Id=I_seq(2)
Ii=I_seq(3)

I0_mod=abs(I0)
I0_ph=angle(I0)
Id_mod=abs(Id)
Id_ph=angle(Id)
Ii_mod=abs(Ii)
Ii_ph=angle(Ii)

%
% Positive sequence circuit
% Thanks to the first circuit in Figure 2, we can see that the positive sequence voltages at points 1 and 2 are: {\bar{E} }_{d,1} ={\bar{E} }_d -{\bar{Z} }_d^S {\bar{I} }_d and {\bar{E} }_{d,2} ={\bar{E} }_d -\left({\bar{Z} }_d^S +{\bar{Z} }_d^l \right){\bar{I} }_d

Ed1=Ed-Zs_d*Id
Ed1_mod=abs(Ed1)
Ed1_ph=angle(Ed1)

Ed2=Ed-(Zs_d+Zl_d)*Id
Ed2_mod=abs(Ed2)
Ed2_ph=angle(Ed2)

%
% Negative sequence circuit
% Thanks to the second circuit in Figure 2, we can see that the inverse sequence voltages at points 1 and 2 are: {\bar{E} }_{i,1} =-{\bar{Z} }_i^S {\bar{I} }_i and {\bar{E} }_{i,2} =-\left({\bar{Z} }_i^S +{\bar{Z} }_i^l \right){\bar{I} }_i

Ei1=-Zs_i*Ii
Ei1_mod=abs(Ei1)
Ei1_ph=angle(Ei1)

Ei2=-(Zs_i+Zl_i)*Ii
Ei2_mod=abs(Ei2)
Ei2_ph=angle(Ei2)

%
% Zero sequence circuit
% Thanks to the second circuit in Figure 2, we can see that the zero sequence voltages at points 1 and 2 are: {\bar{E} }_{0,1} =-{\bar{Z} }_0^S {\bar{I} }_0 and {\bar{E} }_{0,2} =-\left({\bar{Z} }_0^S +{\bar{Z} }_0^l \right){\bar{I} }_0

E01=-Zs_0*I0
E01_mod=abs(E01)
E01_ph=angle(E01)

E02=-(Zs_0+Zl_0)*I0
E02_mod=abs(E02)
E02_ph=angle(E02)

%
% Once we have calculated the values of the voltages using symmetrical components, we can transform the values into the phasor domain.
% For the voltages at point 1: \left\lbrack \begin{array}{c}
% {\bar{E} }_{a,1} \\
% {\bar{E} }_{b,1} \\
% {\bar{E} }_{c,1}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{ccc}
% 1 & 1 & 1\\
% 1 & \alpha^2  & \alpha \\
% 1 & \alpha  & \alpha^2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% {\bar{E} }_{0,1} \\
% {\bar{E} }_{d,1} \\
% {\bar{E} }_{i,1}
% \end{array}\right\rbrack

E1=Ts*[E01; Ed1; Ei1]
Ea1=E1(1)
Ea1_mod=abs(Ea1)
Ea1_ph=angle(Ea1)
Ea1_ph_deg = Ea1_ph*radtodeg;
Eb1=E1(2)
Eb1_mod=abs(Eb1)
Eb1_ph=angle(Eb1)
Eb1_ph_deg = Eb1_ph*radtodeg;
Ec1=E1(3)
Ec1_mod=abs(Ec1)
Ec1_ph=angle(Ec1)
Ec1_ph_deg = Ec1_ph*radtodeg;

spplot([Ea1_mod Ea1_ph*radtodeg],[Eb1_mod Eb1_ph*radtodeg],[Ec1_mod Ec1_ph*radtodeg])
legend('','','','','Ea1','Eb1','Ec1');

% For the voltages at point 2: \left\lbrack \begin{array}{c}
% {\bar{E} }_{a,2} \\
% {\bar{E} }_{b,2} \\
% {\bar{E} }_{c,2}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{ccc}
% 1 & 1 & 1\\
% 1 & \alpha^2  & \alpha \\
% 1 & \alpha  & \alpha^2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% {\bar{E} }_{0,2} \\
% {\bar{E} }_{d,2} \\
% {\bar{E} }_{i,2}
% \end{array}\right\rbrack

E2=Ts*[E02; Ed2; Ei2]
Ea2=E2(1)
Ea2_mod=abs(Ea2)
Ea2_ph=angle(Ea2)
Eb2=E2(2)
Eb2_mod=abs(Eb2)
Eb2_ph=angle(Eb2)
Ec2=E2(3)
Ec2_mod=abs(Ec2)
Ec2_ph=angle(Ec2)

spplot([Ea2_mod Ea2_ph*radtodeg],[Eb2_mod Eb2_ph*radtodeg],[Ec2_mod Ec2_ph*radtodeg])
legend('','','','','Ea2','Eb2','Ec2');